c1.pps=function(data,y,m,M,param){
  data=as.name(data); y=y; m=m
  n=length(m); ybari=y/m
  bias=(1/(n*(n-1)))*((1/n)*(sd(m)/mean(m))^2-cor(m,y)*(sd(m)/mean(m))*(sd(y)/mean(y)))
  plot(m,y,main=data)
  if(param=='mean'){
    est=(1/n)*sum(ybari)
    vhat=(1/(n*(n-1)))*sum((ybari-est)^2)
    bound=2*sqrt(vhat)
    lower=est-bound; upper=est+bound
  } else{
    est=(M/n)*sum(ybari)
    vhat=M^2*((1/(n*(n-1)))*sum((ybari-est)^2))
    bound=2*sqrt(vhat)
    lower=est-bound; upper=est+bound
  } 
  cat("","\n","Results from 1-stage Cluster sample using pps","\n","Data =",data,"\n","n =",n,"N =",N,"\n","Relative bias =",
      bias,"\n","Estimate of",param,"=",est,"\n","Vhat(",param,") =",vhat,"\n","Bound =",bound,"\n",
      "Lower Bound =",lower,"Upper Bound =",upper,"\n","")
  results=list(est=est,data=data,n=n,M=M,vhat=vhat,bound=bound,lower=lower,upper=upper,param=param,bias=bias,y=y,m=m,N=N)
}

# to use the function with its call:
# c1.pps(data,y,m,M,param)

# data: name of dataset, in quotes
# y: numeric vector; the total of all observations in the ith cluster
# m: numeric vector; number of elements in cluster i, i=1,2,...,N
# M: number of elements in population (single number)
# param: c('mean','total')